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SUMMARY 


In the present report a compressible Navier-Stokes solution procedure 
is applied to the flow about an isolated airfoil. Two major problem areas 
have been investigated. The first area is that of developing a coordinate 
system and an initial step in this direction has been taken. An airfoil 
coordinate system obtained from specification of discrete data points has 
been developed and the heat conduction equation has been solved in this 
system. Efforts required to allow the Navier-Stokes equations to be solved 
in this system are discussed. The second problem area is that of obtaining 
flow field solutions . Solutions for the flow about a circialar cylinder and 
an isolated airfoil are presented. In the former case the prediction is 
shown to be in good agreement with data. 



INTRODUCTION 


The prediction of the flow field about isolated airfoils at various 
angles of attack has been a subject of great interest to aerodynamicists 
over the years. The importance of such flow fields is evident in the design 
of aerodynamic machinery. The associated aerodynamic problems vary from a 
whole machine to a detailed analysis of some constituent part. To name a 
few important problems^ it is often necessary to detemnine from the flow 
field the performance of such components as an airplane wing, a control 
surface, a propeller or a helicopter rotor blade. In this context it should 
be clear that a reliable method of predicting aerodynamic flow fields is of 
great value . 

In cases in which the airfoil is at a moderate angle of attack, the 
viscous boundary layer will remain attached over most of the airfoil surface; 
in these situations, lift and moment predictions as well as predictions of the 
detailed pressure distributions along the blade surfaces can be made with 
inviscid flow calculations such as those of Giesing (e.g., Ref. l), Giesing, 
Kalman and Rodden (Ref. 2) or Bauer, Garabedian and Korn (Ref. 3)- If skin 
friction drag for steady flows is desired, a boundary layer procedure such 
as those of Refs. 4 or 5 can be readily performed providing separation does 
not occur. For unsteady flows the incompressible procedures of Nash and 
Patel (Ref. 6) and Briley and McDonald (Ref. 7) are available. If separation 
does occur but the separation region remains small, predictions of the air- 
foil viscous layer can still be made^ as has been demonstrated by Shamroth 
and Kfeskovsky (Ref. 8) and Kreskovsky, Shamroth and Briley (Ref. 9) who 
applied the Briley-McDonald separation bubble calculation method to the 
problem of flow about an airfoil section with varying angle of attack. 

The above procedures are applicable only in the absence of significant 
regions of separated flow. In many cases of practical interest the flow 
about the airfoil contains regions of significant separation? and in these 
cases procedures based upon an outer inviscid flow calculation which ignores 
viscous displacement effects and an inner viscous solution in the immediate 
vicinity of the airfoil are inapplicable. 

An important example of an airfoil flow field which contains significant 
regions of separated flow is the helicopter rotor blade. It is this problem 
which has motivated the present study. Under high speed flight conditions 
the retreating rotor blades are subjected to a diminished dynamic pressure 
and as a result high blade performance requires large lift coefficients to 
be present over the retreating portion of the rotor disc. These large lift 
coefficients are generated by placing the blade at large incidence angles 


3 



relative to the oncoming velocity field. At these large incidences the 
boundary layers developing along the airfoil surface will separate over at 
least a portion of the chord. When separation becomes significant, the blade 
experiences a deterioration in performance and this deterioraticn is termed 
stall. 

Airfoil stall can be divided into two main categories, static stall 
and dynamic stall. Static stall occiars when an airfoil is placed at a large 
incidence angle in a steady stream. This type of stall has been discussed in 
detail by McCullou^ and Gault (Ref. 10). I)ynamic stall occurs when the 
incidence is a function of time. For example, as the helicopter rotor blade 
travels through the rotor disc, it may experience a varying incidence angle. 
Over most of the rotor disc the blade will be unstalled, however, stall may 
occur over a portion of the disc and over this portion the airfoil performance 
will suffer. Eynamic stall differs from static stall in two obvious ways. 
First of all, experimental evidence clearly shows that both the maximum lift 
obtainable and the incidence at which performance first deteriorates are 
greater under dynamic conditions than under static conditions (e.g.. Ref. ll). 
Secondly, although under static conditions lift is uniquely related to 
incidence, the dynamic stall process has a hysteresis loop associated with 
it so that the lift depends upon the incidence history. This historical 
phenomena makes the dynamic stall process much more difficult to predict 
than static stall process. 

Despite its complex nature dsmamic stall has been the subject of several 
recent experimental investigations. For example, the behavior of the leading 
edge separation bubble has been investigated by Velkoff, Blaser and Jones (Ref. 
12) and Isogai (Ref. 13) '-nd the mechanism of dynamic stall on a MCA 0012 air- 
foil has been investigated by McCroskey and Philippe (Ref. lU), McCroskey, 

Carr, and McAlister (Ref. 15 ) end Parker (Ref. l6). The variation of lift and 
moment coefficients through the stall regime has been presented in several 
works such as those of Liiva (Ref. ll) and Landgrebe and Bellinger (Ref. 17). 

Although dynamic stall may occur in a variety of flow situations it 
plays a particularly important role in helicopter rotor performance, partic- 
ularly since a rotor blade may encounter a continuously varying free stream 
velocity and a continuously varying incidence angle as it proceeds around 
the rotor. Obviously the blade performance in terms of lift and moment 
coefficients will depend upon how the blade reacts to its changing environ- 
ment. In addition, blade fatigue stress, blade flutter and aircraft vibra- 
tion will be dependent upon the periodic blade loading and unloading. Thus 
the aerodynamic performance, the rotor structural integrity and aircraft 
vibration characteristics are all significantly affected by possible dynamic 
stall phenomena. 
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To date a variety of approaches have been used to analyze and predict 
dynamic stall. Since the various approaches have been motivated by different 
immediate objectives, the obvious differences in philosophy which have 
appeared are not surprising. A recent review of many of these methods has 
been given by McCroskey in Ref. l8. Ham and Garelick (Ref. I 9 ) and Ham 
(Ref. 20 ) have attempted to model dynamic stall by an inviscid process in 
which vortices are shed from the airfoil leading edge. Although the theory 
has predicted mastlmum lift and moment coefficients during dynamic stall, it 
has not predicted the dynamic stall loop. Furthermore, in this approach the 
time or incidence at which stall begins either must be assumed arbitrarily 
or must be based upon a data correlation. Baudu, Sagner and Bouquet (Ref. 21) 
also have developed an approach based upon a vortex shedding model. 

Other investigators have approached the problem by building semi- 
empirical models of the stall process. These models have been motivated by 
the need for practical prediction techniques and have proven quite useful 
in predicting dynamic stall characteristics. For example, Ericsson and 
Reding have developed a procedure which predicts dynamic stall characteristics 
by combining static airfoil data with semi-empirical models (Refs. 22 and 23). 
Another example of this type of work is that of Lang (Ref. 24) who developed 
a prediction procedure which combines a bubble bursting criteria with an 
inviscid flow analysis. Although these methods serve a pressing practical 
nee<^ their empirical natrore dictates that they be used with caution. Still 
another approach predicts dynamic stall characteristics from data correlations. 
For example. Carta, Commerford and Carlson have related stall characteristics 
to the time-incidence history of the airfoil (Ref. 25). However, since this 
method is based upon data correlation for specific airfoils in specific types 
of motion, it is not clear how this method cotild be extended either to 
alternate airfoil shapes or to other types of motion. 

Althou^ the previously mentioned procedures fill an important immediate 
need, they are limited by semi-empiricism or dependence upon data correla- 
tions. A basic imderstanding of the stall process or the ability to make 
predictions well outside the correlating data base requires more fundamental 
approaches and several such approaches have been developed. For example, 
the analysis of Scruggs, Nash and Singleton (Ref. 26) treats fully turbulent 
boundary layer flow developing under a prescribed pressure gradient obtained 
from the full inviscid equations. The analysis has as its main objective 
an assessment of the effect of unsteady phenomena upon the trailing edge 
separation point. A more comprehensive analysis has been developed by 
Grind, and Reeves (Ref. 27 ) who combined linearized potential flow equations 
with boundary layer equations to predict the flow field behavior. While i 
this procedure has had some success in predicting lift and moment hysteresis 
loops through stall, the model has several shortccmings . First, althou^ the 
Crimi-Reeves approach uses a finite-difference boundary layer calculation in 
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regions of attached flow^ it resorts to an integral procedure in the reversed 
flow region. In addition, the procedure uses a linearized inviscid analysis 
and empirical transition and bubble bursting criteria. In a recent work 
Shamroth and Kreskovshy (Ref. 8) have developed a weak interaction dynamic 
stall analysis which is not limited by integral boundary layer assuicgptions , 
semi-empirical bubble bursting assumptions and semi-empirical transition 
assumptions. In this work Shamroth and Kreskovsky applied the Briley- 
McDonald separated flow finite-difference code (Ref. 7) in conjunction with 
the Giesing potential flow code (Ref. l) to analyze the dynamic stall 
problem. Although the procedure's weak interaction limitation prevented 
the code from calculating hysteresis loops through stall, the analysis did 
predict what appears to be the correct stall mechanism for a MCA 0012 airfoil 
oscillating at a Reynolds number of approximately 0.5 x 107. At the time of 
the Shamroth-Kreskovsky study it was commonly assumed such an airfoil would 
stall as a result of leading edge bubble bursting, however. Ref. 8 predicted 
the bubble to remain attached and stall to originate through abrupt separa- 
tion of the turbulent boundary layer. This stall mechanism has been 
confirmed by the experiments of McCroskey, Carr and McAlister (Ref. 15 ) 
and Parker (Ref. l6), and as pointed out in Ref. 15 such behavior is consis- 
tent with the static data of Gregory and O'Reilly (Ref. 28) and Ridder 
(Ref. 29). 

Despite the fact that both the Crimi-Reeves analysis and the Shamroth- 
Kreskovsky analysis approach the problem from a fundamental basis, both analy- 
ses are still limited by their treatment of the interaction between the inner 
viscous layer and the nominally inviscid flow. The Shamroth-Kreskovsky 
analysis does not allow any interaction and thus cannot be applied in regions 
of significant viscous separation. Although the Crimi-Reeves analysis 
includes an interaction model, this model is based upon a semi -empiric ally 
determined pressure in regions of separated flow. In addition, it should 
be recalled that the Crimi-Reeves inviscid analysis is abased upon linearized 
theory. These considerations indicate the need for new analyses not limited 
by interaction models and an obvious candidate analysis fulfilling this 
requirement would be the full Navier-Stokes equations. 

With the continued improvements in computers and the continued rapid 
advancement in numerical techniques, Navier-Stokes procedures have become an 
increasingly attractive alternative for calculating the flow about an 
airfoil. Much early work in this area concentrated upon predicting flow 
about a circular cylinder and a comprehensive bibliography on this subject 
can be found in the recent paper by Coutanceau and Bouard (Ref. 30). 

It should be noted that those works mentioned in Ref. 30 include solutions 
of both the steady state and time -dependent equations, however, all the 
procedures were limited to incompressible flow. More recently, Navier-Stokes 
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procedures have been applied to more complex shapes including airfoils 
(Refs. 31-32j 34 ). In. Ref. 31 Mehta and Lavan have investigated the flow 
field about an impulsively started airfoil with a finite difference proce- 
dure using stream, function and vorticity as dependent variables. These 
results gave an excellent picture of the starting process. However, the 
Mehta-Lavan procedure is limited to incompressible flow and has been applied 
only in conjunction with a conformal mapping procedTire . In addition, the 
large number of relaxation iterations required per time step may give 
relatively long running times. According to Ref. 32, the run times for 
this method averaged 9-5 minutes of umVAC II 08 time per time step. Despite 
this relatively long rvin time, the excellent results of Ref., 31 represent 
a strong argviment for Navier-Stokes solutions. 

In another analysis, Wu and Sampath (Ref. 32) have applied the Wu- 
Thompson integro-differential formulation (Ref. 33) "to the airfoil problem. 
This calculation also has shed light on the impulsively started airfoil 
problem, however, it is not easily extendable to compressible flow and to 
date it has been formulated and used only for an airfoil which can be mapped 
conformally into a circle. In addition, the procedure requires continuously 
more run time as the rotational portion of the flow field grows. This 
makes the procedure less attractive for predicting steady state solutions 
than for predicting transients. A different approach has been 
taken by Verhoff (Ref. 34) who has applied MacCormack's fully explicit 
method (Ref. 35) to the airfoil problem. Unlike the procedures of Refs. 

31 and 32 , this formulation is compressible, however, being fully explicit 
the procedure requires many time steps to convergence leading to relatively 
long run times. it is likely that this run time problem could be 
alleviated if MacCormack's more recent technique were used (Ref. 36 ). 

Finally, in a recent paper Steger (Ref. 37) has used a viscous analysis 
in conjunction with the coordinate generation procedure of Thompson, Thames 
and Ifestln (Ref. 38 ) to predict flow about an airfoil. The viscous analysis 
is that of Beam and Warming (Ref. 39) which follows Briley and McDonald 
(Ref. 4o) in combining a Taylor expansion linearization with a Douglas-Gunn 
ADI procedure. The major difference in the two approaches lies in the choice 
of dependent variables. The same linearization was used with an alternate 
ADI approach by Lindemuth and Killeen (Ref. 4l). The basic coordinate 
generation procedure (Ref. 38 ) depends upon the solution of an elliptic 
set of partial differential equations. Additional modifications are 
made in one coordinate direction so that mesh points are concentrated 
near the body surface and in the wake region near a specified branch 
cut. The branch cut, however, cannot be technically treated as a 
branch cut in the usual sense. The problem here occurs because only 
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function ‘boundary conditions have "been applied along the cut; the result is 
a cut with tangent discontinuities in the crossing coordinate curves. The 
simultaneous application of both function and derivative conditions would 
certainly remove this problem. However, it would also require the solution 
of a higher order elliptic system: a process which would be costly and 

inefficient. By contrast, the methods developed in this study can effi- 
ciently and correctly model the branch cut with both function and derivative 
specifications . 

The work under the present effort represents the first step in 
obtaining a general efficient computer code capable of predicting flow 
about airfoils of general geometry. The work is an extension of the original 
work of Briley and McDonald (Ref. 4l) for the numerical method and of Eiseman 
(Ref. 42) for the geometry. Although at present only two-dimensional, laminar 
flows are considered, the procedxire is extendable both to three dimensions 
and to turbulent flows . 
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LIST OF SYMBOLS 


Except where dimensions are specified, all quantities in the following 
are nondimens ional; physical velocities are normalized by Up, density by 
p^, pressure by Pj,u2, dynamic viscosity by and time by (L/vl^) where L 
is the reference length. 


a 

A 


b 

B 


c 

ds 



D ,D^ 
m m 

ay 


2 > 

m 



F. . 
xja 


Speed of sound; m/sec 

Constant defined in Eq. ( 13 ); 

Coordinate system parameter, Eq. (15); 

Matrix coefficient, Eq. (25b) 

l^fejor to minor axis ratio of ellipse 

Constant defined in Eq. ( 13 ); 

Coordinate distribution parameter 

Speed of light; m/sec 

Differential element of arc length 

Differential increment of coordinate y^ 

Diameter of circtilar cylinder 

m 

Finite difference operators for coordinate y 
Momentum equation coefficient, Eq. ( 13 ) 

Multidimensional spatial derivative operator vector 

Spatial derivative operator vector associated with coordinate y 

Coordinate tangent vector field, 8x/9y^ 

Momentum equation coefficient, Eq. (12) 

General function, Eq. (20) 

Momentum equation coefficient, Eq. (13) 
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LIST OF SYMBOLS (C0WT*D) 


"id 

^id 


G . , 


h. 




. uo 

■'Q'V 


I 


m 


M 


M 


M 


cyp 


A 

n 


Metric tensor coefficient 
Inverse metric tensor coefficient 
General function, Eq. (20 ) 

Momentum equation coefficient, Eq. (13) 

Metric coefficient, orthogonal coordinates 
General vector function, Eq. (19) 

Momentum equation coefficient, Eq. (13) 
Jacobian 

Momentum equation coefficient, Eq. (13) 

Number of nonlinear equations solved; 
Airfoil chord length 

Reference length, meters 

Momentum equation coefficient, Eq. (13 ) 

Linear operator, Eq. (25c) 

Coordinate distribution parameter, Eq. (16) 

Major axis of ellipse 

Reference Mach number 

Momentxim equation coefficient, Eq. (13) 

Unit normal vector 

Static pressure 
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LIST OF SYMBOLS (COKT’D) 



r 


R 


Re 


s. 

1 


S-| 

-Lmax 


S 

So 

t 

T 

o 
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U. 

1 

A A 




V. 

X 




Momentum equation coefficient^ Eq. (11 ) 

Coordinate parameter, Eq. (l6) 

Coordinate distribution function, Eq. (l6) 

Momentum equation coefficient, Eq. (13) 

Reference Reynolds number, p u L/u 

r r ^ r 

Arc length along coordinate y^ 

Total surface arc length 

General vector function, Eq. (19) 

Momentum equation coefficient, Eq. (13) 

Time; 

Coordinate parameter, Eq. (15) 

Total temperature 

Stress-energy tensor, Eq. (5) 

Streamwise cartesian velocity component 

Physical velocity component along coordinate y^ 

Cartesian unit vectors 

Transverse cartesian velocity component 

Space-time velocity 

Covariant velocity component 

Contravariant velocity component 


11 


LIST OF SYMBOLS (COM"D) 


i 

3tjX 


X 

m 

m 

y.y 

a 


3 

3 

Y 


At 


Ax 


Ay 


m 

m 




0 

y 


p 
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Cartesian position vector 

Arbitrary computational coordinate 

Computational coordinates 

Outer loop position vector, Eq. (l5) 

Spatial differencing parameter, Eq. (29) 

Inner loop position vector, Eq. {ik) 

Time differencing parameter 

Ratio of specific heats, C /C 

P V 


Kronecker delta 


Time increment 

Mesh spacing for coordinate x 

m 

Mesh spacing for coordinate y^ 

Spatial difference operators, Eq. (29) 

Polar coordinate 
rynamic viscosity 

2 

Kinematic viscosity; m /sec 
Density 

General dependent variable 

General dependent variable, Eq. (2U) 

Fourth-derivative dissipation coefficient, Eq. (3^) 



LIST OF SYMBOLS (COWT'D) 


Subscripts 

® Denotes free stream conditions 


r 


Denotes dimensional reference -value 


Superscript 

n Denotes time level n 



ANALYSIS 


The Governing Equations in General Tensor Form 

The role of generalized coordinates . - In the numerical solution of 
fluid dynamic problems there are many advantages to be gained by judicious 
choice of coordinates (eg., Ref. 43). The most obvious advantage is that 
the physical boundaries of a flow region can be represented by coordinate 
surfaces. This removes the need for fractional cells in general; hence, 
the complications and loss of accuracy associated with a boundary inter- 
polation are removed. Another advantage is that a uniform numerical method 
can be used. The solution can then be performed with a fixed number of 
cells in any given direction and with a uniform mesh spacing. The result 
is a simplification of the computer logic; hence, a savings in time for both 
the computer and the programmer results . 

In addition, the coordinate transformation can be constructed to 
contain distributions for physical space mesh points. In this context, the 
unifoim mesh of computational space is simply mapped into a suitably 
distributed mesh in physical space. The resolution of large solution 
gradients is the major objective in the selection of a coordinate mesh 
distribution. A classical example is the resolution of attached boundary 
layers. Another more subtle example is the resolution of large gradients 
in computational coordinates due to regions of high curvatiire on the 
bounding surfaces. When the transformation contains the mesh point distri- 
bution there is no need to construct the apparatus for the discrete approxi- 
mation of derivatives on a nonuniform mesh. This results in a savings in 
both computer logic and storage. As an illustration, consider the case 
where it is desired to automate the difference molecule so that the numerical 
technique can be changed with a few parameters. Changes, in practice, 
usually amount to a selection between forward, backward, or central 
differences. For any given direction, three parameters each for first and 
second derivatives are needed for second order accurate methods. Thus, 
counting 6 parameters for boundary conditions, a total of 12 parameters 
are needed for each ADI direction. This compares favorably with the direct 
approximation of derivatives on a nonuniform mesh where the requirement is 
for 6N parameters on an ADI direction of length N. 

A further advantage of the generalized coordinate approach is that for 
a given problem coordinates can be selected from a large class of cooixiinate 
systems. In the process of sorting through the various possible coordinate 
systems two criteria arise. First, the new coordinates must lead to a real 
simplification; secondly, the coordinates must be easily generated. 



Since bounding surfaces usually become coordinate surfaces the first criter- 
ion is almost always met. The remaining complexity in the first criterion 
is directly measured by consideration of the metric tensor (gij) which is 
obtained from the expression for the fundamental element of arc length 

(ds)^ = Qjj dy' dy* (l) 

Specifically, an increase in the number of nontrivial elements in the 
expression of the metric tensor is accompanied by a corresponding increase 
in the number of terms in the equations of motion. The result is an 
increase in the computational work that is needed after the coordinates 
have been generated along with the necessary metric data. The second 
criterion, \mlike the first, is most often neglected. The unfortunate 
result is that there is often more work involved in making the coordinates 
than in solving the original problem with a less efficient satisfaction of 
the first criterion. In fact, both of the criteria above usmlly are at 
opposite polarities in complexity. The prudent selection of coordinates is 
then a balance between these criteria. 

The criteria for selecting a suitable system of coordinates can be used 
to compare the various classes of coordinate systems and to evaluate the 
relative utility of each. The evaluation will start with conformal trans- 
formations and continually enlarge the class until the most general time- 
dependent coordinates are reached. 

For conformal transformations the metric tensor is simply given by a 
scalar multiple of the identity. That is, gj^j = h(y)6ij where the Kronecker 
S 3 mibol 6j^j vanishes unless i = J in which case it is lonity. From this 
expression it is easy to show that h = where J is the Jacobian of 

the n-dimensional conformal transformation. The simplicity of the metric 
leads to very simple equations of motion at the expense of greatly 
restricting the class of easily obtained transformations. These trans- 
formations are generally obtained by the solution cf partial differential 
equations which may in itself be costly. In addition, the control over 
the mesh distributions is indirect at best. In two dimensions, however, 
conformal transformations have been successfully used on many occasions. 

Here the metric is given by g^j = ( J | 6 j , and the theory of functions of 
one complex variable is a powerful tool which can be used. When the 
boundaries of the flow region can be matched with well-known conformal 
transformations there is nothing that can effectively compete with this 
way of generating coordinates. 

In a number of cases boundaries can be matched through a sequence of 
well-known transformations . However, in most cases of practical importance 
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the boundaries are too coniplicated; and consequently, cannot be siniply 
defined as desired. Thus, one is led to approximate methods. For the 
general airfoil shapes there is the method of Theodorsen (Ref, UU) and along 
similar lines there is the more recent work of Ives (Refs, U5j ^6) which 
employs the Fast Fourier Transform, Both techniques map airfoils to near 
circles through a sequence of well-known maps; and then use a Fourier type 
of approximation. The Schwartz-Christoffel transformation may be used to 
approximate arbitrary shapes with piecewise linear curves. This technique 
works best for simply connected regions where no branch cuts are needed. A 
basic limitation in this method is the poor representation of wall curvature. 
This can be partially resolved by using the Schwartz-Christoffel transfor- 
mation with rounded corners as in Henrici (Ref. ^T). But then there is 
little control over the rounding process. Conformal mappings in higher dimen- 
sions also exist (cf,. Ref. h8) but are generally difficult to construct. 

When conformal mappings become overly difficult to construct, it is 
best to consider the slightly larger class of orthogonal transformations. 

For orthogonal transformations the metric tensor is given by the diagonal 
form gj_j = [hj_(y) Note that, unlike the confonnal transformations, 

the diagonal entries of the metric can be different. The deviation from 
conformality can now easily be measured by an examination of the ratios of 
the functions hj_ as now is demonstrated by an explicit geometric interpre- 
tation of the metric. For a position vector field x, the vector field ej_ = 
(9x)/(8yi) j_s the natural tangent vector field along coordinate curves 
generated by holding the remaining coordinates y^,..., y^"^, y^*^^,..., y^ 
constant. It is often common practice to use the operator notation where 
the position vector field is omitted. By an application of the chain rule, 
the fundamental element of arc length can be expanded as 


(ds)^ = dx.dx 



3x 
d y ' 


ax 

ayi 


dy' dy ^ 


(ej.6j) dy' dyj 


( 2 ) 
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and hence, by linear independence g^j = e^*ej. Now note that the metric is 
orthogonal if and only if the e. and are perpendicular when i n . But 
perpendicularity of e^ and ej at a point is equivalent to the perpendicular 
crossing of the associated coordinate curves at the point in question. 
Consequently the notion of orthogonality in terms of coordinate curves is 
equivalent to the metric expression above. In addition the functions hj^ 
are easily seen to be equal to the lengths of the corresponding natural 
tangent vectors . When the lengths are equal it is clear that the trans- 
fomation is conformal. However, as the ratios of the hj_ deviate from unity, 
the transformation smoothly deviates from conformality. 

With fewer constraints on the metric, the selection of coordinates 
from the class of orthogonal transformations is slightly less restrictive 
than a selection from the class of conformal transformations. The process 
of coordinate generation is usually accomplished by geometric methods. The 
desire is to create families of mutually orthogonal coordinate surfaces. As 
a starting point, one usually begins the process with a given family of 
surfaces that are generated in some way from the boundary of the flow 
region. Families of orthogonal surfaces are then to be constructed to 
complete the specification of coordinates . The first family of s\irfaces 
defines a unique nonnal vector field. This vector field is then extended 
to a smooth field of orthogonal frames which must be integrated to generate 
the orthogonal coordinate surfaces. The condition for integrabillty is 
contained in the Frobenius Theorem (Ref. h9). The computational process 
generally leads to the solution of a system of differential equations. This, 
however, is often a difficult exercise Just to obtain orthogonal coordinates. 

General nonorthogonal coordinates are often preferable to orthogonal 
coordinates since the mesh distributions can be controlled and since the 
coordinates are considerably easier to generate. The construction process 
is usually geometric and generally does not rely on the solution of 
differential equations. Certain methods, however, are not entirely based 
upon the geometry, but upon the solution of a system of elliptic partial 
differential equations (e.g. ,Ref. 38). Such methods are generally of 
ccmparable efficiency with those of the confoimal type. The entirely 
geometric methods usually can be used with distribution fimctions replacing 
independent variables so that mesh point distributions are used directly. 

Such is generally not the case for other methods. The considerable improve- 
ment in flexibility associated with the class of general spatial coordinates 
does come with a small price. Specifically, the metric tensor has 
generally nontrivial off diagonal elements. As with the difference between 
orthogonal and conformal coordinates, the deviation of the general non- 
orthogonal coordinates from orthogonality can be measured directly from the 
metric. That is, the cosine of the angle between distinct coordinate curves 
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is given by the dot product of the associated unit tangent vectors. The 
cosine of the angle between cuives i and j can be written as ; 



gii 


(3) 


Thus when g^j vanishes for distinct i and j the coordinates are orthogonal, 
and when gj^j increases from 0 the coordinates smoothly deviate from 
orthogonality with the deviation given by the arc cosine of Eq. ( 3 ). Local 
regions of nearly orthogonal coordinates can then be constructed within 
a nonorthogonal system. This can be used to some advantage in regions where 
large solution gradients are esipected. For example, borxndary layer 
regions may be treated with nearly orthogonal coordinates which smoothly 
deviate from orthogonality as largely inviscid regions are approached. An 
illustration of this type of construction is given in Ref. 4$. 

When time -dependent problems are considered the general nonorthogonal 
spatial coordinates can be used provided that the boundaries of the flow 
region are rigidly fixed relative to the region. However, if the region 
changes shape as a function of time, then the purely spatial coordinate 
analysis above is no longer valid unless special precautions are taken. 

In terms of the metric the pseudoriemannian metric from special relativity 
is used. Here, the cartesian frame x^, x^, x^ is extended to a Lorentz 
frame by the addition of a time-like coordinate x*^ (Refs. 50-54). The 
fundamental expression for arc length in space-time is 

(ds)^ = dy‘ dyj 

where the summation is now from 0 to 3- 

In this context, the classical equations of motion are obtained from 
the vanishing divergence of the stress-energy tensor and a subsequent 
approximation for slowly moving fluids . This is discussed at more length 
in the following section. Possible applications for such coordinate systems 
include ducts with moving walls such as blood flow problems as well as 
flutter problems where airfoils may be oscillating relative to each other. 

As special cases the classical Eulerian and Lagrangian coordinates as well 
as everything in between are obtained. 

The equations of motion for a viscous fluid . - The equations of motion 
of a viscous fluid are given by a divergence free stress-energy tensor 
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(Refs. 50-54) where v is a space-time velocity, v^ are covariant components, 
are contravariant coiaponents, |a is viscosity, p is pressure, p is density, 
commas denote covariant derivatives, is the inverse metric, and 6^ is the 
Kronecker delta. In tensor analysis, the term contravariant refers to the 
components of a vector in the basis covariant, to the components of the 
vector in the basis defined by the duality condition *^0 = 6^ • Here the 
Latin indices generally vary from 0 to 3 and repeated indices are to be summed 
(Einstein summation convection). The metric is generated from a local Lorentz 
frame and generally depends upon the velocity of light c. Since the veloc- 
ities in classical mechanics are much less than the velocity of light, the 
Navier-Stokes equations can be retrieved as an approximation to the equations 

T |\,=0 ( 6 ) 


P 

when terms of order c“ relative to \mity are removed. The advantage 
inherent in the approximation of the above special relativistic equation 
is that the Navier-Stokes equations are expressed in a manner independent 
of any space-time coordinate system and are given a metric structure induced 
from the relativistic structure. The metric stmicture contains the classical 
Coriolis and centrifugal force effects in a clean and concise manner. That 
is, in addition to spatial changes, the metric contains all of the time- 
dependent variations of the coordinate transformations. The details of 
the approximation were presented in McVittie (Ref. 50 ) for inviscid flows 
and in Walkden (Ref. 52) for viscous flows. All quantities in the following 
equations are nondimensional; physical velocities are normalized by u^., 
density by Pp, pressure by p^Uj,^, dynamic viscosity by |a^, and time by 
(l/uj.) where L is the reference length. The resulting equations for viscous 
flow (Ref. 52 ) are 


— + A 0 (7) 

d t dyP 


for i = 0 and 
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for i = 1, 2, 3 where J is the Jacobian, The equation for i = 0 can be 
identified as the continuity equation; the equations for 1 5. i ^ 3j as the 
respective momentum equations. Here, i = 0 represents the time-like direction 
and Greek indices a, 3j e, to represent space-like directions as they vary only 
from 1 to 3. The energy equation, also presented in Ref. 52, can often be 
replaced with 


= P _ 


A + B 


vi] 


(9) 


under an assumption of constant total temperature, where A = Tq/yM^^ and 

B = -(y-1)/2y. An expansion of the momentum equation leads to a system of the 

form 
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which is suitable for automatic computation. For a time-independent metric, 
this result reduces to the Navier-Stokes equations in a fixed frame 
(Ref. 55). By use of the intermediate quantities 

0 ^ - 9^0 9 °^ ( 11 ) 
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the coefficients of the above system are given by 

K = AJ 


ay 


F|?„ =(gal sf + Bgij hi ) J 

k _ i. iM/S|‘8f + eg.. g‘^) 
2 ay“ \ ^ 'J / 


Ki - _ , ^go/g 

M dp a 

dy 


Rq = A 


AA 


dy 




*^900 . <^gei8 ey3 

■ ^ « nr y 


dy° 


dy' 


'ay ~ Re 


d/! =^J 


I- 8? 8y - gay g^^ — Qy Sq 

o 


~ Re ^..a I 3 ^ dt 2 ^ ^ at 


_ M ‘^gg /3 [ ± _ jL 0*^*0 <^gy 


dy' 


•a I' 


= e“ 

-av '-Qi' ^ ^ 


( 12 ) 


(13) 


U) 


_ -2_ + _A_ /iij _^gae.\ 

3 ay°\f’® at ; ay/9 \ Re at / 

_ aEgy Re ^Qf0 agycu 2 aj g^ysl 

ay^ 2 ay° I ay 3 dy*' J 

[g““ of H. g“^ Q/ - i g*^ 8,“] 


ay' 


21 


The CooMinate System 


An effective coordinate system for the two-dimensional isolated airfoil 
problem can be generated from two loops (or arcs) of data. The two loops first 
are given a parameterization t in a manner which causes both a distribution 
of loopwise points and an alignment of points between the two loops. Here, 
the inner loop is the airfoil; the outer loop, the outer computational 
boxmdary. The coordinate system is then generated by a linear interpolation 
between loops. If Q'(t) denotes an outer loop; P(t), an inner loop; and R(r), 
a distribution function with independent variable 0 < r < 1; then the 
coordimate transformation is given by 


X = /§'+ R{a-p) 


( 14 ) 


where the cartesian location x is 3 when R = 0 and 5^ when R = 1. The 
metric data for this coordinate system is presented in Ref. k2 where the loops 
are constructed for a cascade of airfoils rather than for an isolated airfoil. 
Although the coordinate transfomation, Eq. (l4), is 'Representative of a 
large class of transformations ^ it is instructive to specialize our 
discussion at least initially for purposes of illustration. Thus, in most 
of the discussion, it will*'be assimed that the airfoil has no cusps and is 
convex; minor modifications can he employed for cusps and regions of 
concavity. With convexity, an orthogonal coordinate system can be generated- 
if each interpolation between inner dnd outer loop points is over a line 
normal to the airfoil surface. This can be accomplished by generating the 
outer computational boixndary u as a uniform dilation along airfoil normal 
lines^ That is, if n(t) is the outward unit normal vector field to the air- 
foil P(t), then the outer computational boundary is given by 

a(t) =^(t) + Afi(t) (15) 

for some fixed distance A. The parameter t, as a loopwise distribution, 
can easily by specified as the arc length of some intermediate loop "^(t) + 
Bn(t) where 0 ^ B ^ A. The parameter t can be viewed as a label for 
corresponding inner and outer loop points. This is analogous to inner and 
outer points being displayed as sequences each with an index t. In taking 
the analogy one step further, the pseudoradial curves correspond to the 
lines joining inner and outer loop points with the same index. As B is 
adjusted from 0 to A the distribution varies from airfoil arc length towards* 
outer loop arc length. Under a uniform discretization of t, this will 
cause points to cluster around the airfoil regions of higher curvature 
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(notal)ly at leading and/or trailing edges) as A is approached. To complete 
the specificaiton of the coordinate transformation, Eq. (lU), a definition of 
R(r) is needed. The function R(r) is used to specify a distribution of 
points along the coordinate curves normal to the airfoil. Commonly, it is 
desired to resolve a boundary layer near the airfoil surface. For this 
purpose, the distribution 


R (r) = mr + (I- m) 


tgnh P(l-r) 
tanh D 


(16) 


is most effective (Ref. h2) , Here, a slope m is selected to determine a 
suitable uniform boundary layer mesh corresponding to the line R = mr; the 
level of adherence to this line is a result of the selected damping factor D. 
As D is increased from 0, the distribution will depart from global uniformity 
and more closely follow the given line before increasing to the end point 
R(i) = 1. By construction, the only inflection point occurs at R(i) = 1. 

Due to the inflection point condition, the distribution leaves an almost 
uniformly distributed boundary layer region with a continuously expanding 
grid as r approaches 1. 


Once the coordinate transformation is established, the numerical 
solution for a desired flow field can be attempted. The mesh in the 
computational domain is given by uniform discretizations of both r and t. 

The array of mesh points here is a rectangular foimation which is consider- 
ably simpler than the corresponding array in physical space. The advantages 
of this formulation are clear. 


In certain cases the coordinate system can be constructed directly from 
an analytically specified airfoil surface. For example, consider an airfoil 
in the shape of an ellipse. The analytic specification is given by 


iS(0) = 


M (cos 0 U| + sin 9 Ug) 
■v/! +( — |) sin 6 


(17) 


for 0 ^ 6 ^ 2tt and where u^ and U 2 are the standard cartesian unit vectors . 
The major axis is denoted by M; the ratio of major to minor axes, by b. If 
the outward expansion along the unit nomnals is sufficiently large (as is 
often the case), then the outer loop is nearly circular. As such, the 
distribution t determined by the outer loop arc length is, for all practical 
puiposes, proportional to the angles formed by the outward \onit normals 
and the x - axis . With a little algebra, it can be shown that 
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( 18 ) 


e = ,a„-' 

which upon substitution leads to the desired parametric form of the airfoil 

l(t). 


Generally, the coordinate construction must be done numerically. The 
process is broken down into a sequence of operations; namely, the determina- 
tion of airfoil normals, the selection of B, the construction of inner and 
■outer loops , and finally the generation of coordinates from R and the two 
loops. In more detail these operations are tabulated as follows: 

(a) Generate the airfoil surface 3(s) where s is an arc length 
parameteri zat ion . 

(b) Construct the unit vector field n(s). 

(c) Choose B such that o < B< A 

(d) Construct the intermediate loop Y-n(s) = 3(s) + B n(s) 

(e) Store the paired points (^(s), y^(s)) 

(f) Reparameterize by its arc length t to obtain yg('t) 

(g) From the pairing in (e), reparameterize the airfoil surface to 
obtain ^(t) 

(h) Generate the outer boundary from Eq. (l5) 

(i) Choose the number, N, of loopwise mesh points 

( j ) If t varies from o to some value T, then set At = T/(N-1) 

(k) For each t = nAt , evaluate inner and outer loops as n = o,..., N-1 

(i 1 Choose the parameters m and D which control the desired resolution 
for the boundary layer region 

(m) Choose the number, M, of mesh points in the pseudoradial direction 

Cn) For each r = m/(M-l), evaluate the radial distribution (Eq, (l6)) 
for m = o, 1, . . . , M-1 

(o) Generate an NxM coordinate grid from Eq. (lU), step (k) and step (n) 
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If the airfoil data is not smooth, then a least squares curve fit will he 
needed to determine the normals (Ref. 56). Further curve fits will be needed 
to construct inner and outer loops. If, however, B is not too large, then 
inner and outer loops can be constructed directly from a fit to the inter- 
mediate curve corresponding to B. Since the intermediate curve has an arc 
length parameterization it is usually easier to obtain an accurate fit . In 
any case, it is only necessary to fit one curve; the others can be obtained 
as expansions or contractions along the normal lines. 

The method of coordinate construction presented here can also be 
adapted to the specific high Reynolds number problem considered by Steger 
(Ref. 37). In that problem the mesh was generated to resolve a flow field 
about an isolated airfoil with a long but narrow wake region. For this case, 
the leading edge region on the airfoil is considered to be bounded by the 
points above and below the airfoil where the outward unit normal vector first 
becomes vertical. The unit vector field ii(t) is now chosen to coincide with 
the outward unit normal vector in the leading region, to point vertically 
outward in the remaining airfoil regions, and in continuation to point 
vertically away from a horizontal branch cut in the wake region. The rdsult 
is two arcs of data rather than two loops of data. As before, the distribu- 
tion t is determined by the arc length of an intermediate arc. The general 
construction and analysis of two loop (or arc) coordinate systems is given 
in Ref. 4$. There the analysis includes coordinate stretches which are 
applicable to the wake region where an expanding mesh may be desired. An 
illustration of this coordinate system is given in Fig. 1. 

The MIKT Procedure 

One of the major obstacles to the routine nmerical solution of the 
miilti -dimensional compressible Navier-Stokes equations is the large amount 
of computer time generally required, and consequently, efficient computa- 
tional methods are highly desirable. Most methods for solving the compres- 
sible Navier-Stokes equations have been based on explicit difference 
schemes for the vinsteady form of the governing equations and are subject to 
one or more stability restrictions on the size of the time step relative to 
the spatial mesh size . These stability limits usually correspond to the 
well known Courant.-Friedrichs-Lewy (CFL^ condition and, in some schemes, 
to an additional stability condition arising from viscous teims. In one 
dimension, the CFL condition is At ^ Ax/(|u| + a), and the viscous stability 
condition is At s (Ax)^/2v, where At is the time step. Ax is the mesh size. 
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u is the velocity, a is the speed of sound, and v is kinematic viscosity. 
These stability restrictions can lovrer computational efficiency by imposing 
a smaller time step than would otherwise be desirable. Thus, a key dis- 
advantage of conditionally stable methods is that the maximum time step is 
fixed by the spatial mesh size rather than the physical time dependence or 
the desired temporal accuracy. In contrast to most explicit methods, 
implicit methods tend to be stable for large time steps and hence offer 
the prospect of substantial increases in computational efficiency, provided 
of course that large time steps are acceptable for the physical problem of 
interest and that the computational effort per time step is competitive 
with that of explicit methods . In an effort to exploit these potentially 
favorable stability properties, an efficient Implicit method based on 
alternating-direction differencing techniques was developed by Briley and 
McDonald (Ref. 4o). 

This method was subsequently designated the Multidimensional Implicit 
Nonlinear Time -dependent (MINT) solution procedure. The MINT method was 
further developed and applied to both laminar and turbulent duct flows by 
Briley, McDonald and Gibeling (Ref. 57). Subsequently, a three-dimensional 
compressible Navier-Stokes combustor flow analysis employing the MINT 
procedure was developed by Gibeling, McDonald and Briley (Ref. 58) and this 
procedure was then employed by Levy, et al., (Ref. 59) to determine the 
feasibility for computing a turbulent shock-wave boundary layer interaction 
with an implicit method. 

Outline of method . - The MINT procedure has been previously described 
in Refs, 4o and 57j however, the description will be repeated here for 
completeness. The method can be briefly outlined as follows: the governing 

equations are replaced by an implicit time difference approximation, option- 
ally a backward difference or Crank-Nicolson scheme. Terms involving 
nonlinearities at the implicit time level are linearized by Taylor expansion 
about the solution at the known time level, and spatial difference approxi- 
mations are introduced. The result is a system of m\alti dimensional coupled 
(but linear) difference equations for the dependent variables at the 
unknown or Implicit time level. To solve these difference equations, the 
Douglas-Gunn (Ref. 60) procedure for generating alternating-direction 
implicit (ADI) schemes as perturbations of fundamental implicit difference 
schemes is introduced. This technique leads to systems of coupled linear 
difference equations having narrow block-banded matrix stmctures which 
can be solved efficiently by standard block- elimination methods. 

The method centers around the use of a formal linearization technique 
adapted for the integration of initial-value problems. The linearization 
technique, which requires an implicit solution procedure, permits the 
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solution of coupled nonlinear equations in one space dimension (to the 
requisite degree of accToracy) by a one-step noniterative scheme. Since 
no iteration is required to compute the solution for a single time step, 
and since only moderate effort is required for solution of the implicit 
difference equations, the method is computationally efficient; this 
efficiency is retained for multidimensional problems by using ADI techniques. 
The method is also economical in terms of computer storage, in its present 
form requiring only two time-levels of storage for each dependent variable. 
Furthermoare, the ADI technique reduces multidimensional problems to sequences 
of calculations which are one-dimensional in the sense that easily-solved 
narrow block-banded matrices associated with one-dimensional rows of grid 
points are produced. Consequently, only these one-dimensional problems 
reqviire rapid-access storage at any given stage of the solution proced-ure, 
and the remaining flow variables can be saved on auxiliary storage devices 
if desired. 

Althoiogh present attention is focused on the compressible Navier-Stokes 
equations, the numerical method employed is quite general and is formally 
derived for systems of governing equations which have the following form: 

where 0 is a column vector containing i dependent variables, H and S are 
column vector functions of 0, and .2) is a column vector whose elements are 
spatial differential operators which may be multidimensional. The generality 
of Eq. (19) allows the method to be developed concisely and permits various 
extensions and modifications (e.g., noncartesian coordinate systems, 
turbulence models) to be made more or less routinely. It should be empha- 
sized, however, that the Jacobian 8H/80 must usually be nonsingiilar if the 
ADI techniques as applied to Eq. (19) are to be valid. A necessary condition 
is that each dependent variable appear in one or more of the governing 
equations as a time derivative. An exception would occur if for instance, 
a variable having no time derivative also appeared in only one equation, 
so that this equation co\ild be decoupled from the remaining equations and 
solved a posteriori by an alternate method. As a consequence, the present 
method is not directly applicable to the incompressible Navier-Stokes 
equations except in one dimension, where ADI techniques are \jtnnecessaiy. 

For example, the velocity-press\are form of the incompressible equations has 
no time derivative of pressiare, whereas the vorticity-stream- function form 
has no time derivative of stream function. For computing steady solutions, 
however, the addition of suitable "artificial" time derivatives to the 
incompressible equations, as was done in Chorin's (Ref. 61) artificial 
compressibility method, woiald permit the application of the present method. 
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Alternatively, a low Mach number solution of the present velocity-density 
formulation of the compressible equations can be computed. 

Linearization technique . - A number of techniques have been used for 
implicit solution of the following first-order nonlinear scalar equation in 
one dependent variable 0(x,t): 

a<^/dt =F(<^) aG(<^)/ax (20) 


Special cases of Eq. (20) include the conservation form if F(0) = 1, and 
quasi-linear form if G(0) = 0. Previous implicit methods for Eq. (20) 
which employ nonlinear difference equations and also methods based on two- 
step predictor-corrector schemes are discussed by Ames (Ref. 62, p. 82) 
and von Rosenburg (Ref. 63, p. 56). One such method is to difference 
nonlinear terms directly at the implicit time level to obtain nonlinear 
implicit difference equations; these are then solved iteratively by a 
procedvire such as Newton's method. Although otherwise attractive, there 
may be difficulty with convergence in the iterative solution of the non- 


linear difference equations, and some efficiency is sacrificed by the need 
for iteration. An implicit predictor-corrector technique has been devised 
by Douglas and Jones (Ref. 64) which is applicable to the quasilinear case 
(G = 0) of Eq. (20). The first step of their procedure is to linearize the 
equation by eyaluating the non- linear coefficient as F(0^) and to predict 
values of 0^‘'’2 using either the backward difference or the Crank-Nicolson ^ 
scheme. Values for 0^'''^ are then computed in a similar manner using F(0^"’'^) 
and the Crank-Nicolson scheme. Gourlay and Morris (Ref. 65) have also 


proposed implicit predictor-corrector techniques which can be applied to^ 

Eq. (20). In the conservative case (F=l), their technique is to define G(0) 
by the relation G(0) = 0G(0) when such a definition exists, and to evaluate 
^(0n+l) using values for 0^'*‘^ computed by an explicit predictor scheme. With 
$ thereby known at the implicit time level, the equation can be treated as 
linear and corrected values of are computed by the Crank-Nicolson scheme. 


A technique is described here for deriving linear implicit difference 
approximations for nonlinear differential equations. The technique is based 
on an expansion of nonlinear implicit terms about the solution at the known 
time level, t^, and leads to a one-step two-level scheme which, being linear 
in unknown (implicit) quantities, can be solved efficiently without iteration, 
This idea was applied by Richtmyer and Morton (Ref. 66, p. 203) to a scalar 
nonlinear diffusion equation. Here, the technique is developed for problems 
governed by i nonlinear equations in i dependent variables which are func- 
tions of time and space coordinates. Although the present effort concen- 
trates upon two spatial dimensions and time, the technique will be described 
for the three-dimensional, unsteady equations . 
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The solution domain is discretized by grid points having equal spacings 
in the computational coordinates, Ay^, Ay2 and Ay3 in the y^, y^ and y3 
directions, respectively, and an arbitrary time step. At. The subscripts 
i, j, k and superscript n are grid point indices associated with y^, y , y^ 
and t, respectively, and thus ip2,jfk 'denotes 0(y^, yj, y^, • I't is assumed 

that the solution is known at the n level, t^, and is desired at the (n+l) 
level, At the risk of an occasional ambiguity, one or more of the 

subscripts is frequently omitted, so that 0^ is equivalent to j i,. 


The linearized difference approximation is derived from the following 
implicit time-difference replacement of Eq. (19) ; 




( 21 ) 


where, for example, s The form of .25 and the spatial differenc- 

ing are as yet unspecified. A parameter (O ^ P ^ l) has been introduced so 
as to permit a variable centering of the scheme in time. Equation (21) 
produces a backward difference formulation for P = 1 and a C rank- Ni cols on 
fomulation for p = ■§■. 


The linearization is performed by a two-step process of esspansion about 
the known time level t’^ and subsequent approximation of the quantity 
(p0/9t)^At, which arises from chain rule differentiation, by (0’^'*'^-0*^) . 

The result is 


hH+i _ + 0(At)^ 

sn+'=s" + (dS/<3<^)^ (<^" + '-<^")+0(At)^ 


(22a) 

(22b) 

(22c) 


The matrices 8H/80 and 8S/80 are standard Jacobians whose elements are 
defined, for example, by (811/80)^^. = 811^/80^. The operator elements of 
the nBtrix 8.2>/90 are similarly ordered, i.e., (8.2)/90)qj. = 8.25 however, 

the intended meaning of the operator elements requires some clarification. 

For the qth row, the operation (8.2>q/90)^ _ 0^) is understood to mean 

that |9/8t.2)g^[0(x,y, z,t)]} ^ At is computed and that all occurrences of 
(80^/8t)^ arising from chain rule differentiation are replaced by 
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(0j+l _ 0^)/At. 

After linearization as in Eqs. (22), Eq. (21) 'becomes the following 
linear implicit time-differenced scheme; 

Although is linearized to second order in Eq. (22a), the division by 

At in Eq. ( 2 l) introduces an error term of order At. A technique for 
maintaining formal second-order accuracy in the presence of nonlinear time 
derivatives is discussed by McDonald and Briley (Ref. 4o), hovrever, a three- 
level scheme results. Second-order temporal accuracy can also be obtained 
(for 3 = ■!■) t)y a change in dependent variable to 0 s H(0), provided this 
is convenient, since the nonlinear time derivative is then eliminated. 

The temporal accuracy is Independent of the spatial accuracy. 

On examination, it can be seen that Eq. (23) is linear in the quantity 
( 0 .n+l _ 0 n) that all other quantities are either known or evaluated at 
the n level. Computationally, it is convenient to solve Eq. (23) for 
( 0 n+l _ 0 tt) rather than 0^"''^. This both simplifies Eq. (23) and reduces 
roundoff errors, since it is presumably better to compute a small 0(At) 
change in a 0(l) quantity than the quantity Itself. To simplify the 
notation, a new dependent variable ii defined by 

0 = < 0 _ 0 n (24) 


is introduced, and thus = 0 ^'''^ - 0 ^, and i//^ = 0. It is also convenient 

to rewrite Eq. (23) in the following simplified form: 

(A+At/)'/'"'^' = At [.Si (<^")+s"] 

where the following symbols have been introduced to simplify the notation: 

As (3H/a<^)" -/5At(dS/dc/))^ (2Jb) 


i’s (25c) 
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It is noted that l{ii) is a linear transformation and thus I{0) = 0. 
Furthermore, if 2>(0; is linear, then ^(0) = -32>(0). 

Spatial differencing of Eq, (25a) is accomplished simply by replacing 
derivative operators such as 3/9y“, 9^/9y“9y™ by corresponding finite 
difference operators, Djjj, Henceforth, it is assumed that ^ and J. 

have been discretized in this manner, unless otherwise noted. 

Before proceeding, some general observations seem appropriate. The 
foregoing linearization technique assumes only Taylor esjpandability, an 
assumption already implicit in the use of a finite difference method. 

The governing equations and boundary conditions are addressed directly as 
a system of coupled nonlinear equations which collectively determine the 
solution. The approach thus seems more natural than that of making ^ hoc 
linearization and decoupling approximations, as is often done in appl3d.ng 
Iniplicit schemes to coupled and/or nonlinear partial differential equations. 
With the present approach, it is not necessary to associate each governing 
equation and boundary condition with a particular dependent variable and 
then to identify various "nonlinear coefficients" and "coupling terms" 
which must then be treated by lagging, predictor-corrector techniques, 
or iteration. The Taylor expansion procedure is analogous to that used in 
the generalized Newton-Raphson or quasi-linearization methods for iterative 
solution of nonlinear systems by expansion about a known current guess of 
the solution (e.g.. Bellman & Kalaba, Ref. 67). However, the concept of 
expanding about the previous time level apparently had not been employed to 
produce a noniterative implicit time -dependent scheme for coupled equations, 
wherein nonlinear terms are approximated to a level of accuracy commensurate 
with that of the time differencing. The linearization technique also permits 
the implicit treatment of coupled nonlinear boundary conditions, such as 
stagnation pressure and enthalpy at subsonic inlet bouindarles, and in prac- 
tice, this latter feature was found to be crucial to the stability of the 
overall method (Refs. 4o and 57). 

Application of alternating-direction techniques . - Solution of Eq. (25a) 
is accomplished by application of an alternating-direction implicit (ADl) 
technique for parabolic -hyperbolic equations. The original ADI method was 
introduced by Peaceman and Rachford (Ref. 68) and Douglas (Ref. 69); 
however, the alternating-direction concept has since been expanded and 
generalized. A discussion of various alternating-direction techniques is 
given by Mitchell (Ref. 70) and Yanenko (Ref. 7l). 

The present technique is simply an application of the very general 
procedure developed by Douglas and Gunn (Ref. 60 ) for generating ADI schemes 
as perturbations of fundamental implicit difference schemes such as the 
backward-difference or Crank-Nicolson schemes . 
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For the present, it vri.ll be assumed that .2)(0) contains derivatives of 
first and second order vfith respect to y^, and y^, but no mixed deriva- 
tives. In this case, can be split into three operators, 2>2» 

associated vd.th the y^ and y^ coordinates and each having the functional 
form .2)^ = ^(0, 9/9y™, 9^/9y®^9y®). Equation (25a) then becomes 

[A+At(i', +^2 + ^3 At[(^, +2>2 + ^3)‘^"+S"] (26) 

Recalling that 2(}lP') = 0, the Douglas-G'unn representation of Eq. (26) can 
be vrritten as the follovri.ng three-step solution procedure: 

(A+At^, )'/'* = At[( +^2 + 2>3)<^"+s"] ( 27 a) 


( A + At i’g )<//** = aV'* (27b) 


(A + At =AV'** (27c) 

where if and 0 are intermediate solutions. It will he shown subsequently 
that each of Eqs. (27) can be written in narrow block-banded matrix form 
and solved by efficient block-elimination methods. If 0 and 0 are 
eliminated, Eqs. (27) become 

(A + At/, )A-'(A + At 72)A-'(A + At ^ = At [( <#>"+s"] (^8) 


If the multiplication on the left-hand side of Eq. (28) is performed, it 
becomes apparent that Eq. (28) approximates Eq. (26) to order (At)2. 
Although the stability of Eqs. (27) has not been established in circum- 
stances sufficiently general to encompass the Navler-Stokes equations, it 
is often suggested (e.g., Rlchtmyer & Morton, Ref. 66, p. 215) that the 
scheme is stable and accurate under conditions more general than those 
for which rigorous proofs are available. This latter notion was adopted 
here as a working hypothesis supported by favorable results obtained in 
actual computations (e.g.. Refs. 40, 57-59)* 

A major attraction of the Douglas-Gunn scheme is that the intermediate 
solutions and 0** are consistent approximations to Furthermore, 

for steady solutions, 0^ = 0 =0 = 0“ independent of At. Thus, 

physical bo\mdary conditions for 0^"^1 can be used in the intermediate 
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steps -without a serious loss in accuracy and -with no loss for steady 
solutions. In this respect, the Douglas-Gxmn scheme appears to have an 
advantage over locally one-dimensional (LOD) or "splitting" schemes, and 
other schemes whose intermediate steps do not satisfy the consistency 
condition. The lack of consistency in the intermediate steps complicates 
the treatment of ho-undary conditions and, according to Yanenko (Ref. 71j 
P« 33 )j does not permit the use of asymptotically large time steps. It is 
not clear that this advantage of the Douglas-G-unn scheme would always 
outweigh other benefits which might be derived from an alternative scheme. 
However, since the ADI scheme can be viewed as an approximate technique 
for solving the fundamental difference scheme, Eq. (25a), alternate techni- 
ques can readily be used -within the present formulation. 

It is worth noting that the operator can be split into any number 
of components which need not be associated -with a particular coordinate 
direction. As pointed out by Douglas and Gunn (Ref. 6 o), the criterion 
for identifying sub-operators is that the associated matrices be "easily 
solved" (i.e., narrow-banded). Thus, mixed derivatives can be treated 
implicitly within the ADI framework, although this would increase the 
mmiber of intermediate steps and thereby complicate the solution procedure. 
Finally, only minor changes are introduced if, in the foregoing development 
of the numerical method, H, .2>, and S are functions of the spatial 
coordinates and time, as well as 0 . 

Solution of the implicit difference equations . Since each of Eqs. (27) 
is implicit in only one coordinate direction, the solution procedure can 
be discussed -with reference to a one-dimensional problem. For simplicity, 
it is sufficient to consider Eq. (27a) -with ^ 2 ? -^3 ~ Consider the 
following three -point difference formulas : 

D^<#>s[a A_ + (i-a)A^.]<#*/Ax^= (a<#>/dx^)j + 0 [Ax^+(a-l/2) Ax^]( 29 a) 


= (d^«^/ax^i+ 0 (Ax^ ( 29 b) 


for a typical computational coordinate Xj^. Here, A_ = 0^ - 0i_i> ^+ = 
0j^^2. ~ ^±> parameter a has been introduced (O s ot s l) so as -fco 

permit continuous variation from backward to forward differences. The 
standard central difference formula is recovered for or = ^ and was \ised 
for all solutions reported here. 
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As an example, suppose that the qth component of the vector operator ^ 
has the form 

,2 


m 


= -a7-6„(+) +F/,(*) G„(*) (30) 

^ ^ .m 


where F and G are column vector functions having the same (hut an arbitrary 
number) of components; F^ denotes the transpose of F. The form of Eq. (30) 
permits governing equations having any number of first and second spatial 


A, '^■^'"<1 


Then, 
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(^n+i 


[f^ ^ 
Liq dx„ 
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(31) 


It is now possible to describe the solution procedure for Eq. (27a) 
for the one-dimensional case with given by Eq. (30) and difference 
formulas given by Eq. (29). Because of the spatial difference operators, 

Djjj and Eq. (27a) contains ijj* ;//. , and consequently, the 

system of linear equations generated by writing Eq. (27a) at successive 
grid points (Xjjj)j_ can be written in block- tridiagonal form (simple tridia- 
gonal for scalar equations, Ji = l) . The block-tridlagonal matrix structure 
emerges from rewriting Eq. (27a) as 






(32) 


where a, b, c are square matrices and d is a column vector, each containing 
only n- level quantities. When applied at successive grid points, Eq. (32) 
generates a block-tridiagonal system of equations for tjy which, after 
appropriate treatment of boundary conditions, can be solved efficiently 
using standard block- elimination methods as discussed by Isaacson and 
Keller (Ref. 72, p. 58)* The solution procedure for Eqs. (27b, c) is 
analogous to that just described for Eq. (27a). It is worth noting that 
the spatial difference parameter of can be varied with i or even term by 
term. For example, an "upwind difference" formula can be obtained if a is 
chosen as 1 or -1 depending on the sign of the elements of F^; however, the 
formal accuracy of the method would then be reduced to first order. 
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Computing requirements , - Various block-elimination algorithms can be 
devised for solution of equations with block-banded matrix stiuctures (cf., 
Isaacson & Keller, Ref. 72). Such algorithms can be derived using variants 
of Gaussian elimination for a banded matrix, but with the square submatrix 
elements of the banded matrix processed using matrix algebra. Thus, opera- 
tions involving matrix subelements are not assumed to commute, and division 
by a matrix subelement is accomplished by computing the inverse and multi- 
plying. Following this procedure, McDonald and Briley (Ref. 4o) have 
developed an algorithm for block-tridiagonal systems arising from the second- 
order difference formulas, Eq. (29)- The algorithm requires only one 
inverse per grid point. A standard operation count (scalar multiplications 
and divisions) has been performed for systems with L x L block elements and 
N diagonal block elements, i.e., L coupled equations along IT grid points. 

The block-tridiagonal scheme requires (3N-2) (l 3 + L^) operations, the same 
as the matrix factorization scheme of Isaacson and Keller (Ref. 72). 

Ass\jming there are N grid points in each coordinate direction, the total 
number of operations for a single time step is obtained from the operation 
count for solution of one block-banded system by multiplying by 2IT and 31^ 
for two and three dimensions, respectively. 

For comparison, it is noted that in the case of the Navi er- Stokes 
equations, merely evaluating the right-hand side of Eq. (27a), which would 
be a minimum requirement for a one-step explicit scheme, requires 302 N 
operations for a 3-point difference formula. 

In view of the many factors involved, it is difficult to evaluate 
precisely or with any generality the overall computational efficiency of 
the MINT method relative to various other methods. However, the foregoing 
operational counts show that the effort expended to solve the implicit 
difference equations by block-elimination is not excessive compared with 
that necessary simply to evaluate the differenced Navier-Stokes equations, 
let alone the various other bookkeeping tasks present in most large-scale 
computer programs for fluid dynamics problems. In the solutions presented 
here, the solution of the block-tridiagonal systems using double precision 
arithmetic required only about twenty percent of the total computer time 
per time step. 

Artificial dissipation . - In computing solutions for high Reynolds 
number flows, it is often necessaiy to add a fomn of artificial viscosity 
or dissipation. Artificial dissipation in some form is often useful in 
practical calculations to stabilize the overall method when function boundary 
conditions are applied, when coarse mesh spacing is used, or in the presence 
of discontinuities. The need for artificial dissipation arises in certain 
instances when centered spatial difference approximations are used for first 
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derivative terms. The use of artificial dissipation is thus a matter of 
spatial differencing technique, and it is commonly employed in either 
explicit or inherent form, and in both explicit and implicit difference 
schemes . 

One possible dissipation term in common use is based on an observation 
(e.g., Roache, Ref. 73> p. l62) that for a linear model problem representing 
a one-dimensional balance of convection and diffusion terms, solutions 
obtained using central differences for the convection term are well behaved 
provided the mesh Reynolds number ReAxj^ = (■Ujul Ax^ Re is s 2, but that 
qualitative inaccuracies (associated with boundary conditions) occur for 
^®Axra ^ This suggests the use of an artificial viscosity term of the 
form ejj^D^0, where 



to ensure that the local effective mesh Reynolds number is no greater than 
two. The simplified analysis presented here was extended for the generalized 
tensor equations presented previously and the resulting dissipation teims 
were added to the continuity and momentum equations . 

A second type of artificial damping which is a fourth-order dissipation 
term has been suggested by Beam and Warming (Ref. 39) "to damp small wave- 
length disturbances . In the present formulation an explicit fourth-order 
damping term was added directly to the fundamental difference scheme, 

Eq. (25a), as follows; 


(A + At./) 


s"]+ i 

m=i 


d <p 
8 (3x ^ 


(34) 


Note that the dissipative coefficients 0 )^^^ implicitly contain the time step 
At; however, since the dissipation term is treated explicitly (to retain the 
block tridiagonal matrix structure), a linear von Neumann stability analysis 
shows that the coefficients must be in the range of 0 s ^ 1 for 
stability (Ref. 39) • This condition actually Implies a limit on the 
maximimi permissible time step which may be taken, however, in the present 
formulation this was not a restriction. The advantage of the fourth- 
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derivative dissipation added in Eq. (3^) over the more conventional 
artificial viscosity formulation, Eq. (33) > is that the formal accuracy of 
the method is not altered, -whereas Eq. (33) reduces the formal accuracy to 
first order when R eA-ir^ > 2. Strictly speaking, the overall method is 
second-order accurate since -♦ 0 as the mesh is. refined. It should 

be remembered, however, that such asymptotic truncation error estimates are 
meaningful only for sufficiently small mesh size; whereas in practical 
calculations of complex flows mesh resolution capabilities are almost 
al-ways strained. 
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RESULTS 


Solution of the Heat Conduction Equation 

Solutions of the heat conduction equation were the first test cases 
of the MIKE code applied in its general tensor form. In vector notation 
the equation is given by 


v2 (/> = di>/d\ 


(35) 


which can he expressed in tensor form as 



(36) represents the i^^ general coordinate and g. . are the 
o “^ficients of the metric tensor. It should he noted that Eq. (36) is not 
a trivial test case in as much as the equation contains hoth metric 
coefficients and their derivatives. 

Two cases were considered. In hoth cases the solution domain was 
the region interior to the closed curve formed hy an elliptical arc, the 
x-axis and an arc generated hy points with an equal normal distance from the 
elliptical arc. An illustration is shown in Fig. 2. In hoth cases the 
coordinate system was generated hy the two arcs consisting of the elliptical 
arc and the outer arc . The elliptical arc was given as raw geometric data 
and then fit with a parametric curve developed in a piecewise fashion. The 
fitting process used least squares so that eventually measured data could he 
used. The parameterization for the outer arc was taken as its arc length; 
for the elliptical arc, it was imposed from the outer arc under the condition 
that lines joining points of like parameter values are normal to each arc. 
Given this parameterization, the coordinate transformation is given in 
Eq. (l4). 

In the first test case initial conditions were specified such that all 
boundary points were equal to zero and all interior points were equal to 
unity. The zero values on the houndaiy were held as hoimdaiy conditions. 

The solution then decayed to zero throughout the domain as the equation 
was marched in time. In the second case zero derivative conditions were 
applied at the inner and outer curved boundaries. The segments of the x-axis 
were set at zero and unity, respectively. Initially all interior points 
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were set to unity and the solution was then allowed to develop in time. 

At large times the solution satisfied the condition = 0 as it should. 
These solutions of the heat conduction equation do not represent as 
stringent a test of the coordinate system as would a solution of the Navier- 
Stokes equations. However, the heat conduction equation is not trivial 
and the solutions obtained demonstrate, the potential of the current 
coordinate system. A feature of this system is that it does not require 
analytically specified cu27ves as boundaries. Solutions of the Navier-Stokes 
equations also were attempted in this same coordinate system which was 
generated from raw geometric data; however, in this case successful solutions 
were not obtained. 

The source of the difficulty in the Havier-Stokes case was the appear- 
ance of second derivatives of the metric coefficients. In contrast the heat 
conduction equation only contained first derivatives of the metric coeffi- 
cients. "When the boundary was expressed as a non- analytic ally specified 
curve, first derivatives of the metric coefficients were sufficiently 
smooth to obtain reasonable solutions, however, second derivatives were not. 

Present efforts now are aimed at resolving this problem; two approaches 
are being taken. In the first approach the equations are being reformulated 
in a manner which does not require second derivatives of the metric coeffi- 
cients. In the second approach the curve fitting routines are being refined 
further so as to produce sufficiently smooth derivatives to allow a Navier- 
Stokes solution. 


Flow About a Circular Cylinder 

The first Navier-Stokes solution calculated under the present effort 
was the flow about a circular cylinder at Reynolds number (based upon 
diameter) equal to forty. This case was chosen since it represents a 
relatively simple geometry and since both experimental data and other 
solutions were available for comparison, A comprehensive list of references 
on both experimental work and numerical calculations is given by Ref. 30, 

The dependent variables used for this calculation were the density and 
the contravariant components of the velocity, v^ and v^. The calculation 
was initiated by specifying the inviscid solution throughout the flow field 
and then imposing a zero velocity condition on wall boundary points. The 
velocity was set equal to zero and the inviscid transverse momentimi equation 
(8p/8n = O) was solved at the cylinder surface. The latter bo\mdary condi- 
tion is required to determine a value of density at the wall grid point and, 
in the present case of constant total temperature, this condition is 
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equivalent to 8p/9n = 0. It should he noted that the condition of constant 
total temperature and zero slip velocity implies 9T/9n = 0 at the vail. At 
the outer boundary, velocity and density were set to their inviscid values 
over the upstream three quarters of the outer boundary (O G <3 tt/U). First 
derivatives of the physical velocity components vere set to zero and the 
pressure was set to its inviscid value over the remainder of the far field 
boundary. The far field boundary vas taken to be fifteen diameters from the 
cylinder center (A = 29)* In this computation the only artificial viscosity 
term added vas that based on the mesh Reynolds number criterion, Eq. (33). 

The present predictions vere obtained with a 35 x 35 mesh (including boundary 
points), and the coordinate system employed is shown in Figs. (3a,b). A 
special case of Eq. (l6) is the Roberts boundary layer transformation 
(Ref. T8) -which is obtained by setting m = 0, With a damping factor D = 2.7, 
this transformation was used to concentrate grid points both near the vail and 
near the front and rear stagnation points as shown in Fig. 3. The loopvise 
distribution vas computed by a rigid translation of Eq. (l6) with m = 0. 

The result was the distribution 0(0) = tanh (2eD/ir -D)/tanh D where D = 1.5- 
The computation time for the present nonorthogonal form of the governing 
equations, Eqs . (7, 9-13), is approximately 5-9 x 10"^ CPU minutes per grid 
point per time step on the UNIVAC 1110, and approximately 80 time steps vere 
required to obtain the steady state solution. 

The present prediction of the surface pressure along with predictions 
of Son and Hanratty (Ref. T^) and Kavaguti (Ref. 75) are shown in Fig. h. 

The Kavaguti pressure prediction vas given relative to the predicted rear 
stagnation point pressure and in this figure the Kavaguti rear stagnation 
point pressure vas arbitrarily set at the Son and Hanratty value. The 
predictions of Refs. 7^ and 75 both vere obtained from solutions of the 
incompressible Navier-Stokes equations. The present solution is obtained 
from the compressible equations at a Mach number equal to 0.2. As shown in 
Fig. h all three solutions are in reasonable agreement with each other. 

The discrepancies may be due to the different numerical methods or the effect 
of compressibility. Another possible source of discrepancy in the leading 
edge region is the application of the inviscid transverse momentum equation 
as a wall boundary condition in the present formulation. An improved wall 
boundary condition is probably required for proper representation of both 
the leading edge stagnation region and separated flow regions. It is 
believed that a one-sided difference representation of the normal momentum 
equation at the wall is a more realistic boundary condition, since 
this boundary condition allows a normal pressure gradient consistent 
with the momentum balance in the near wall region. The centerline velocity 
prediction is compared to the data of Coutanceau and Bouard (Ref. 30) 
and the predictions of Kawaguti (Ref. 75)9 Apelt (Ref. 76 ) and Nieuwstadt 
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and Keller (Ref. 77) in Fig. 5- Again the agreement between data and all 
analyses is good. Predicted streamline locations are shown in Fig. 6 and 
streamwise velocity profiles at several azimuthal locations are shown in 
Fig. 7. A comparison of Figs. 5 and 6 shows a discrepancy in the wake 
length, however, this discrepancy represents less than one half the radial 
grid spacing at the near wake closure point. The values of the wake length 
given by other authors (Refs. 30, 75^ 76, 77) are higher than the present 
prediction, however, the present results were obtained with a relatively 
coarse mesh. Also, use of the contravariant velocity components, which was 
subsequently abandoned in the airfoil computations, is believed to increase 
the truncation error in the numerical predictions. 


Flow About a Joukowski Airfoil 

The second Navier-Stokes solution calculated was the flow about an 
11 percent thick Joukowski airfoil at zero angle of attack and a Reynolds 
number of eighty (8o) based on the airfoil chord length. The case repre- 
sents a much more severe test of the MINT tensor code because of the signif- 
icant variation in the metric tensor coefficients throughout the computa- 
tional field. The coordinate systems for the Joukowski airfoil computations 
presented herein were obtained using the analytic Joukowski transformation 
in conjunction with Eq. (l6) (for m = 0 and D = 3-3) in the pseudoradial 
coordinate direction. No redistribution of mesh points in the aximuthal 
(loopwise) direction was required for the calculations presented. The 
Joiikowski airfoil results at Re = 80 were obtained using a 4l x 22 mesh 
(including boundary points) with the far field boundary taken to be approxi- 
mately fifteen chord lengths from the airfoil. The coordinate system 
employed is shown in Figs. (8a,b). 

The dependent variables employed in this calculation were the density 
and the physical velocity components, u^ and U 2 * The calculation was 
initiated by specifying the inviscid solution throughout the flow field and 
then imposing a zero velocity condition on wall boundary points. The 
velocity was set equal to zero and the inviscid transverse momentum equation 
(0p/9n = O) was solved at the airfoil surface. At the outer boundary, 
velocity and density were set to their inviscid values over the upstream 
three quarters of the outer boTindary. First derivatives of the velocity 
components were set to zero and the press-ore was set to its inviscid value 
over the remainder of the far field boundaiy. In this computation, both 
the fourth-order damping, (Eq. 34), and the artificial viscosity term based 
on the mesh Reynolds mmiber criterion, (Eq. 33) ^ were added to the governing 
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equations, "because it was found that the conventional dissipation term alone 
(Eq. 33) did not prevent spatial oscillations from occurring. The relatively 
coarse grid employed in' the pseudo- radial direction (22 points) led to some 
difficulties at the far field boundaries where the velocity components were 
specified. Specifically, the computed velocities near the boundary did not 
approach the inviscid values as uniformly as desired. 

The predicted airfoil surface pressure distribution is shown in Fig. 9 
and computed streamline locations are shown in Fig. 10. The computed stream- 
wise and transverse velocity profiles are shown in Figs. 11 and 12, respect- 
ively, at the azimuthal locations defined by the nondimens ional surface arc 
length coordinate, s^/s^^ (see Fig, 8). The computed centerline velocity 
profile downstream of the trailing edge is shown in Fig. 13. These predic- 
tions generally exhibit qualitatively reasonable behavior, except in the 
vicinity of the airfoil leading edge where the pressure seems to be some- 
what higher than expected. This behavior might be attributable to the 
application wall boiuidary condition 9p/9n = 0, which is physically incorrect 
in the leading edge region. 

Finally, the Navier-Stokes calculation for the flow about an 11 percent 
thick Joukowski airfoil at zero angle of attack and a Reynolds number of 1000 
was initiated to demonstrate the high Reynolds number capability of the MIRT 
code. The Joukowski airfoil results at Re^^ = 1000 were obtained using a 
4l X 30 mesh (including boundary points) with the far field boundary taken 
to be approximately ten chord lengths from the airfoil. The coordinate 
system employed is similar to that shown in Fig. 8, but with a higher 
concentration of grid points near the airfoil surface. This calculation 
was begun by applying approximate boundary layer corrections to the inviscid 
velocity field. Again the velocity was set equal to zero and the inviscid 
transverse momentum equation (9p/9n = O) was solved at the airfoil surface. 
This calculation was not run to convergence because of time and funding 
constraints. The surface pressure prediction is shown in Fig. l4, the 
streamline pattern is shown in Fig. 15 and the computed streamwlse velocity 
profiles at several azimuthal locations are shown in Fig. l6 for this case. 
Again, these resvilts appear to be qualitatively reasonable although the flow 
over the airfoil surface had not yet separated at the time of these plots. 



CONCLUSIOUS 


Two major areas have been investigated under the present effort. The 
first area concentrates upon the development of a coordinate system for 
calculating the flow about an isolated airfoil. The second area concentrates 
upon the flow field calculation itself. 

A fundamental requirement for any airfoil flow calculation procedure is 
the generation of an efficient and accurate system of coordinates . The 
efficiency is needed if eventual applications to three-dimensional and/or 
time-dependent geometries are contemplated; the accuracy is needed to 
insure the fiondamental integrity of the overall algorithm and to control its 
error growth. 

Under the present effort the process of coordinate generation for the 
entire class of airfoil geometries has progressed in an orderly fashion 
leading to a solution of the heat conduction equation in a coordinate 
system generated from a discrete specification of boundary points. The 
major problem preventing a solution of the Navier-Stokes equations with 
this coordinate generation procediu*e appears to be a lack of smoothness in 
the higher derivatives of the metric coefficients . Under investigation are 
both an alternate form of the Navier-Stokes equations requiring less smooth- 
ness in metric data and a new coordinate generating technique with stronger 
smoothness properties. 

The second portion of the current effort has successfully solved the 
Navier-Stokes equations about isolated bodies. Predictions for low Reynolds 
number flow about a cylinder are in good agreement with both experimental 
data and other numerical predictions. Calculations of the flow field about 
a Joukowski airfoil at Re^^^ = 80 are qualitatively reasonable; at Re^ = 1000 
no adverse high Reynolds number effects are observed. 

Several problems which have arisen in the course of the present effort 
require further work. These include a thorough investigation of both wall 
and far field boundary conditions. As mentioned previously a one-sided 
difference representation of the normal momentum equation at the wall 
shoiild be considered as a boundary condition for determining the wall value 
of density. In the far fields solution of the Euler equations over a portion 
of the outer computational boundary may alleviate the numerical difficulties 
encoiintered in that region of the flow field, without resorting to mesh 
refinement in the inviscid portion of the flow field. 
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Figure 1.— A coordinate system for flows with a narrow wake. 






Figure 30* - Circular cylinder coordinate system, Re^ = 40; near field (radial locations 1 through 17). 
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Figure 3b* - Circular cylinder coordinate system, Re^ = 40; far field (radial locations 17 through 35)* 
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Figure 5. — Centerline velocity distribution behind circular cylinder, ReQ=40, 
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Figure 7. - Strearavri.se velocity profiles about circular cylinder, Re^ = 4o. 
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Figure 8b. - Joukowski airfoil coordinate system. 




= 80; far field (pseudo-radial locations 12 through 22). 









Nondimensional arc length along coordinate lines 



Figure 12. - Transverse velocity profiles about Joukowski airfoil, Re^ = 80. 
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